Total forces in the diffusion Monte Carlo method with nonlocal pseudopotentials 
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We report exact expressions for atomic forces in the diffusion Monte Carlo (DMC) method when 
using nonlocal pseudopotentials. We present approximate schemes for estimating these expressions 
in both mixed and pure DMC calculations, including the pseudopotential Pulay term and the Pulay 
nodal term. Harmonic vibrational frequencies and equilibrium bond lengths are derived from the 
DMC forces and compared with those obtained from DMC potential energy curves. Results for four 
small molecules show that the equilibrium bond lengths obtained from our best force and energy 
calculations differ by less than 0.002 A. 
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I. INTRODUCTION 

The diffusion Monte Carlo (DMC) method is the most 
accurate approach available for calculating total ground- 
state energies of large many-electron systems^. Energy 
derivatives, and in particular atomic forces, are of great 
importance as they may be used to relax atomic struc- 
tures and perform molecular dynamics simulations. It 
has, however, proven difficult to develop accurate and ef- 
ficient ways of calculating atomic forces within the DMC 
method. 

The DMC method involves using the imaginary-time 
Schrodinger equation to project away excited states of 
a many-electron wavefunction so that the ground-state 
wavefunction $ remains. The fermionic symmetry is 
maintained by the fixed-node approximation^ which con- 
strains the nodal surface of $ (the hypersurface on which 
$ is zero) to equal that of an antisymmetric trial wave- 
function ^Pt- The standard DMC algorithm generates 
the "mixed" probability distribution ^J't'I' which can be 
used to calculate unbiased expectation values (apart from 
the fixed-node error) of operators that commute with the 
Hamiltonian, H . If, however, an operator does not com- 
mute with H , the "pure" probability distribution 4>$ is 
required which can be generated within the DMC method 
using, for example, the future-walking algorithm^. 

Within the Born-Oppenheimcr approximation the 
atomic positions arc treated as parameters rather than 
dynamical variables, and the total atomic force is de- 
fined as the negative energy gradient with respect to the 
atomic position. In the mixed and pure DMC methods, 
the total force estimators are constructed by taking the 
gradient of the (mixed and pure) DMC energy. As in 
other electronic structure methods, these estimators con- 
sist of contributions from the Hellmann-Feynman (HFT) 
forced and from additional Pulay terms^ which must 
be included to obtain unbiased estimates of the forces. 
Although the estimators for the HFT force have simi- 
lar forms in the mixed and pure DMC methods, there 
are significant differences between the Pulay terms. The 
mixed DMC Pulay term involves the derivative of the 
unknown DMC wavefunction $, which cannot be cal- 



culated straightforwardly. Reynolds et al^^ suggested 
using the derivative of ^t instead and obtained a prac- 
tical, although approximate, scheme for estimating the 
mixed DMC Pulay term. The pure DMC Pulay term in- 
volves an integral over the nodal surfacoiS, which cannot 
be calculated in a standard DMC calculation. A prac- 
tical but approximate scheme for estimating this nodal 
surface term has recently been developed^. 

Although mixed DMC total forces have been investi- 
gated in studies of small molecule a^^'^^i^^ , the pure DMC 
nodal surface ter m^^i^^ has not been calculated for real 
systems. In a recent study of a model system without 
electron-electron interactions, however, the nodal term 
was found to be important^. Also, no previous study 
has directly compared the performance of the mixed and 
pure DMC forces for real systems. 

It is well-known that the HFT estimator has an infinite 
variance when the bare Coulomb potential is used to de- 
scribe the electron-nucleus interaction. Different routes 
have been proposed to address this problem. Assaraf 
et a/iiiii^ii^ added a term to the HFT force which has a 
zero mean value but greatly reduces the variance of the 
estimator. Chiesa et alr^ developed a method to filter 
out the part of the electron density that gives rise to the 
infinite variance. Using soft pseudopotentials also elimi- 
nates the infinite variance problemi^, and this method is 
used in the current work. 

Pseudopotentials not only resolve the infinite variance 
issue when calculating forces, they also remove the chem- 
ically inert core electrons and their rapid spatial varia- 
tions from the problem. This greatly reduces the com- 
putational cost of a DMC calculation which scales with 
the atomic number, Z, as ^5.5^ i7,i8 However, the use of 
pseudopotentials introduces additional Pulay-like terms 
in the HFT force estimator— which have been neglected 
in previous calculations. In this work we have included 
these pseudopotential Pulay terms. 

We investigate the accuracy of the mixed and pure 
DMC force estimates for the H2, SiH, GeH and SiH4 
molecules. Bond lengths and harmonic vibrational fre- 
quencies are obtained from the forces and are com- 
pared with those obtained from DMC energy calcula- 



tions. These results are used to both evaluate the im- 
portance of the different Pulay terms and to compare 
the performance of the mixed and pure DMC force esti- 
mators. 

This paper is organized as follows. In Sec. 2 we give ex- 
act and approximate expressions for the forces under dif- 
ferent pseudopotential localization schemes. In Sec. 3 we 
describe our DMC calculations and in Sec. 4 we present 
and discuss the molecular bond lengths and vibrational 
frequencies obtained. We draw our conclusions in Sec. 5. 



II. FORCES IN THE DIFFUSION MONTE 
CARLO METHOD 

We write the valence Hamiltonian for a many-electron 
system as 



H — Hioc 



W, 



(1) 



where -ffioc consists of the kinetic energy, the Coulomb 
interaction between the electrons and the local pseudopo- 
tential, and W is the nonlocal pseudopotential opera- 
tor. Two different pseudopotential localization approx- 
imation (PLA) schemes have been introduced to eval- 
uate the nonlocal action of W on the DMC wavefunc- 
tion, <&. In these schemes H is replaced by an effective 
Hamiltoniani^i2£, 



and upon which both the nodal surface (via ^t) and the 
DMC wavefunction $ depend. Taking the derivative of 
the DMC energy with respect to A gives 



dEp 
dX 



^H'^ + ^(h - Eb) *' 



dV 



$*dV 



$' [H - Ed] * 



dV 



(5) 



^^dV 



for both the mixed and pure DMC methods. We use 
the notation ct' = -^ where a can be a function or an 
operator. The first term in Eq. ([5]) is the HFT force^ 
and the others are Pulay terms^. 



Mixed DMC forces 



^tot 
mix' 



obtained by setting ^P = 'I't in Eq. ^. After some 
rearrangements, we arrive at 



ptot 
niix 



PHFT 

mix 



mix 



mix 



mix' 



(6) 



with 



Ha = Hi 



loc 



iy*i 



Hb = gioc + ^ +W- 

Wt 



(2) 
(3) 



The nonlocal pseudopotential operator W'^ corresponds 
to all positive matrix elements (r'| W~^ |r'), and W~ cor- 
responds to all negative matrix elements^, where r is 
the 3iV-dimensional position vector for the N electron 
system and N is the total number of electrons. Follow- 
ing Ref.— , these two approximations are referred to as 
the fuU-PLA (FPLA) and semi-PLA (SPLA) when using 
Ha or Hb, respectively. The corresponding fixed-node 
DMC ground-state wavefunctions are denoted by $^ and 

The DMC energy can be written in the form 



Er 



^H^dV 



$*dy 



(4) 



which includes the mixed DMC (^ = ^t) and pure DMC 
(\l/ = <i>) estimates of the energy. In all later expressions, 
$ stands for either <i>^ or ^b , and H for either Ha or Hb ■ 
Although the mixed and pure estimates of E^ in Eq. (g]) 
are equivalent for a given localization approximation, E-£, 
may differ under the two localization schemes. 

We now consider a general parameter A, e.g., a nu- 
clear coordinate, which is used to vary the Hamiltonian, 



7.HFT 



F' 



F^ 



F'^ 






$*tV^' dV^ 



Rn 



^^rdV 



E ^> 



Rr 



/3- 



f^ifS^a) 



Ra-R/3|3 



^Vj/n 









dV 



(7) 



(8) 



(f>Vl/r, 



^' \H - Ed] *t 



$ 



*n 



dV 



(9) 



<j)V]/r, 



f&Vj/rpdT/ 



H-Ed]^'t 



*T 



dV 



(10) 



^•^TdV 



-^mix i^ the mixed DMC HFT force and the other expres- 
sions are Pulay terms. The HFT force in Eq. ([7]) contains 
two contributions from the pseudopotential, one from its 
local part l^oc and one from its nonlocal part W, and a 
third contribution from the nucleus-nucleus interaction. 



In this nucleus-nucleus interaction term, Rq represents 
the 3-diniensional position vector of the ath nucleus, and 
Za is the associated charge. The three Pulay terms in 
Eqs. ([5))- (|10p are identified as follows: Fj^i^ results from 
the PLA and is therefore called the pseudopotential Pu- 
lay term, F^^^ is the volume term, and F^;^ is called 
the mixed DMC nodal term since it can be written as an 
integral over the nodal surfaco^. Note that all terms in 
Eqs. (|7|)- (fT0l) take the same form under both localization 
schemes, the only difference is the distribution [^'^^^ or 
^t*&b) used to evaluate the expectation values. A sim- 
ple way to understand this is to note that iJ always acts 
on the trial wavefunction ^I^t and Ha^t = Hb'^t- 

In mixed DMC simulations, it is straightforward to 
evaluate the contributions to the force except for the 
volume term -F^jx, because it depends on the derivative 
of the DMC wavefunction, $'. Since it is unclear how 
to evaluate $' in mixed DMC calculations, we use the 
Reynolds' approximatior^^ 









(11) 



which is exact on the nodal surface (see Eqs. (4) and (16) 
of Ref.— ) but introduces an error of first order in (^t — 'f) 
away from the nodal surface. 



K 



pure 



IVr^l^'dS* 



$$dV" 



(15) 



F^^^^^^ is the pure DMC HFT force, F^^^.^ is the pure 
DMC pseudopotential Pulay term, and the pure DMC 
nodal term F^^^^ is an integral over the nodal surface F 
(defined by ^^t = 0). Where terms appear in braces, the 
upper one refers to the FPLA and the lower to the SPLA. 
The form of the nodal term in Eq. P^ is independent 
of the localization scheme. The nodal term involves the 
gradient Vr$ evaluated at the nodal surface P. Rcfji 
shows that this gradient is defined as its limit when ap- 
proaching the nodal surface from within a nodal pocket 
(where $ is nonzero). The derivation of the nodal term 
from Eq. (O can be found in RefsF^ii^. 

Although the HFT force F^JJ under the FPLA can be 
calculated in the pure DMC method, it is not straight- 
forward to evaluate the action of the nonlocal operator 
{W~y on the unknown DMC wavefunction $b in i^j^^J 
under the SPLA scheme. Therefore, the following local- 
ization approximation. 



(W-)'$B (T^-)'*T 



<^B 



*T 



(16) 



B. Pure DMC forces 



The total force in the pure DMC method, F^°l^, is 
obtained by setting ^ = $ in Eq. ([5]). After some ma- 
nipulations, we obtain 

(12) 



ptot _ pHFT I pP _|_ pN 
pure pure "*" pure "*" pure' 



with 



F, 



HFT 






$A*AdV^ 






jw-y^B 



dV 



$$K,dF 



<^<^dV 



$B<f>sdl/ 



^. E ^.#^^#^ (13) 



/3(a#/3) 



$A$^ 






W^t\ *t 



*n 



Vj/r, 



dV 



-F = < 

pure * 



$B$ 



B'i'B 



^A^AdV 



V]/r, 



\ *T 



V]/r, 



, 1 (14) 
dV 



^B^BdV 



is used in the evaluation of K^^J under the SPLA scheme 
which introduces an error of first order in (\['x — *&). 

Another complication arises with the pure DMC nodal 
term F^^^^.^ in Eq. (jTS)) because it involves the evaluation 
of an integral over the nodal surface. It is unclear how to 
evaluate such an integral in a standard DMC simulation. 
The following relationship suggested in Refii, 



pure mix 



^[(^T-*) 



(17) 



allows the approximate evaluation of F^^^.^ as twice its 
mixed DMC counterpart while introducing an error of 
second order in (^t — *&)• F^^ix ni^Y be evaluated in a 
standard DMC simulation using the volume integral rep- 
resentation of Eq. pO|) . Equation pT]) is an application 
of the standard extrapolation teclmiquoi, as in this case 
the variational estimate of the nodal term is zero^. 

It is worth stating that the pseudopotential Pulay 
terms in both mixed and pure DMC simulations vanish 
when ^T equals $, which follows by inspection. Also, 
the mixed and pure DMC nodal terms vanish when the 
nodal surface of ^t is exact. The proof of this statement 
is analogous to the one presented in Appendix C of Refpi. 



C. Total versus partial derivatives 

Since the atomic force equals the negative total deriva- 
tive of the DMC energy with respect to a nucleus posi- 
tion A, all previous expressions involve total derivatives, 
in particular ^^. Calculating the total derivatives in- 
volves knowledge of how the wavefunction changes with 



A. In the variational Monte Carlo (VMC) method, all 
parameters {q} that describe the wavefunction \E'x can 
in principle be uniquely specified, e.g., by minimizing the 
variational energy. The specification of all {c;} parame- 
ters does not, however, uniquely determine the deriva- 
tive of ^T with respect to A. In standard quantum 
chemistry methods, the derivatives of the Ci with re- 
spect to A arc typically obtained by second-order per- 
turbation theory2ii2^. Unfortunately, such an approach 
is not straightforward in VMC and DMC methods. 

In this work, we follow a different route and approx- 
imate all total derivatives by their partial derivatives, 
which introduces an error of first order in (^t — ^)- We 
expect, however, this approximation to be rather accu- 
rate for the following reason: taking the total derivative 
of the DMC energy with respect to A gives 



dEu 
dX 



dEu 
dX 



EdEi) dci 
~d^~dX' 



(18) 



where the Ci are the parameters in ^t and the Hamil- 
tonian. The sum in Eq. (J18p stems from the implicit 
dependencies of the parameters Ci on A. This sum is ne- 
glected when all total derivatives are replaced with par- 
tial derivatives in all previous force expressions. Since the 
DMC energy is approximately minimized with respect to 
the Ci, we expect that the parameters c; have little ef- 
fect on the DMC energy, i.e., ^o. jg small. Therefore, 
neglecting the sum in Eq. (jlSp . or cquivalently replacing 
all total derivatives with partials in our previous expres- 
sions, is expected to be a good approximation. 



III. DETAILS OF QMC CALCULATIONS 

We use a trial wavefunction of the standard single- 
determinant Slater- Jastrowi form. The orbitals forming 
the Slater-determinants are obtained from Hartree-Fock 
calculations using the GAMESS-US^ code with atomic- 
centered Gaussian basis sets. For all molecules, we use a 
basis set of sextuple-C quality (without / and g functions 
but with four additional diffuse p and d functions). To 
study the influence of the basis set we also use a smaller 
Gaussian basis set for the SiH and GeH molecules with 
only five s and two p- functions so that the trial wavefunc- 
tions ^T and the nodal surface are less accurate. We refer 
to these two basis sets as large and small. 

Table U shows that, when using the small instead of 
the large basis set, the DMC total energies increase 
by 1.01(3) mHa and 2.15(3) mHa for the SiH and GeH 
molecules, respectively. 

We use Jastrow factors consisting of electron-electron, 
electron- nucleus, and clcctron-electron-nucleus terms, 
which are expanded in natural power series^ii. The wave- 
function for H2 has 87 variable parameters, while those 
for the other molecules have 157. The parameters in 
the Jastrow factors are optimized by first minimizing the 
variance of the local energy^^, and subsequently minimiz- 
ing the energ y^^i^^ while keeping the cutoff parameters of 



the Jastrow factor fixcd^l. We use Dirac-Fock averaged 
relativistic effective pseudopotentials22i2^ which can be 
obtained online^. The future-walking method^ is used 
to sample the pure estimates and all DMC calculations 
are performed using the CASINO code^ version 2.1. 

We use the analytic expressions derived in Ref.— for 
evaluating the HFT force. The Pulay terms may also 
be evaluated using analytic expressions, but to make the 
code more easily adaptable to other trial wavefunction 
forms we use a finite-difference approach. This intro- 
duces an error which is linear in the infinitesimal nuclear 
displacement, A. We find that A w 10~^ A minimizes 
the resulting error in the Pulay terms which is about 
seven orders of magnitude smaller than the estimated 
values of the total forces. 

DMC calculations suffer from systematic errors arising 
from the short-time approximation to the Green's func- 
tion, which we have carefully investigated. We find that 
the forces calculated with DMC timesteps of 0.01 and 
0.005 a. u. agree with forces obtained from extrapolations 
to zero timestep within one standard error of less than 
0.0005 a. u. Figure [1] shows such an investigation for the 
SiH molecule using the small basis set where the forces 
acting on the Si and H atoms are plotted as a function 
of timestep. We therefore use a timestep of 0.005 a.u. for 
all our DMC calculations, to avoid the necessity for re- 
peated extrapolation to zero timestep. It is worth noting 
that, for the systems studied here, the timestep errors in 
the HFT and Pulay forces tend to cancel one another. 
This can be seen for the SiH molecule in Figure [1] when 
comparing the solid lines (total forces) with the dashed 
lines (HFT forces). 

For the future-walking pure DMC estimates to be ac- 
curate, an infinite future-walking projection time is in 
principle required. However, we found a projection time 
of 10 a.u. to be sufficient in our calculations as no sig- 
nificant changes in the estimates were found when using 
longer projection times. For example. Figure [5] shows the 
convergence of the pure HFT forces with respect to the 
projection time for the SiH molecule and both basis sets. 

To determine the equilibrium bond lengths, we calcu- 
late the forces at 0%, ±1.5%, ±3%, and ±4.5% around 
the experimental bond lengths. We then fit the deriva- 
tive of the Morse potentials^ with three free parameters 
to our force data and locate the zero-force geometry and 
the harmonic vibrational frequency. For all molecules, we 
compare results derived from the Morse potential with 
those obtained from quadratic and cubic fitting forms. 
We find that the influence of the different fitting forms 
on the derived bond lengths and frequencies is much less 
than one standard error except for a few cases when the 
forces are obtained from the mixed DMC total force es- 
timator. There, the statistical error is much reduced 
and the influence of the fitting form on the derived bond 
lengths and frequencies is sometimes as large as two stan- 
dard errors. Table |lTl for example, presents results ob- 
tained from the different fitting forms as calculated for 
the SiH molecule with both basis sets. It is not obvious 



BL (A) Basis gHF(Ha) gvMc(Ha) EpMcCRa.) Eg 

Ha 0.740 large -1.13367 -1.17399(1) -1.17452(1) 98.7% 

SiH 1.520 large -4.26235 -4.36967(4) -4.37719(2) 93.4% 

GeH 1.600 large -4.24392 -4.34377(1) -4.35143(2) 92.9% 

SiH4 1.480 large -6.08924 -6.27247(3) -6.27927(7) 96.4% 

SiH 1.520 small -4.24689 -4.36563(6) -4.37611(2) 91.9% 

GeH 1.600 small -4.23275 -4.34156(2) -4.34928(2) 93.4% 

TABLE I: Hartree-Fock (HF), VMC and DMC energies (Ha) for four molecules. The first column specifies the bond length 
(BL) used, the second column gives the basis set, and the last column states the percentage of the DMC correlation energy 
retrieved within the VMC method, Ec, as a measure of the accuracy of the Jastrow factor as defined in Refii^. The error bars 
in Ec are smaller than 0.1 %. 



Form Basis a 



P(3/2) 
P(4/3) 
Morse 



small 1.5242(7) 
small 1.5238(7) 
small 1.5242(6) 



1.5138(1) 
1.5138(1) 
1.5141(1) 



1.5259(8) 
1.5261(7) 
1.5259(6) 



2000(18) 
1983(55) 
1992(12) 



2096(4) 
2095(5) 
2084(2) 



2050(14) 
2018(28) 
2045(11) 



P(3/2) 
P(4/3) 
Morse 



large 1.5195(8) 
large 1.5194(8) 
large 1.5195(7) 



1.5105(2) 
1.5107(2) 
1.5107(1) 



1.5175(11) 
1.5177(11) 
1.5173(9) 



2069(18) 
2104(56) 
2049(13) 



2089(4) 
2078(6) 
2080(2) 



2061(18) 
2046(38) 
2052(11) 



Expt. 



1.520 



1.520 



1.520 



2042 



2042 2042 



TABLE II: Bond lengths, a (A), and harmonic vibrational frequencies, to (cm~^), for the SiH molecule derived from different 
functional forms previously fitted to the DMC energies and forces. Three fitting forms are used as indicated in the first column, 
P(3/2) indicates that a cubic polynomial is fitted to the energies and a quadratic one to the forces, and similarly for P(4/3), 
and Morse indicates that a Morse potential is fitted to the energy and its derivative to the forces. The results are compared 



with the experimental values taken from Ref.- 



which fitting form is best suited to describe our fitted 
data. However, the Morse potential gives the smallest 
statistical error bars for all derived bond lengths and fre- 
quencies, which suggests that the Morse potential may 
indeed give the best description of our data. We also ob- 
tain bond lengths and frequencies from a Morse potential 
fitted to the calculated DMC energies. 

The statistical error bars in the calculated bond 
lengths and frequencies are determined using a statistical 
method. For convenience, the DMC forces and energies 
evaluated at the seven different bond lengths are called 
a set of data. For each set of data, statistical noise is 
added where the noise is obtained from a Gaussian dis- 
tribution with a variance specified by one standard error 
of the DMC force and energy estimates. 10^ different 
sets of data are generated in this manner and a Morse 
potential is fitted to these sets. The bond lengths and 
frequencies are obtained for all sets, they are averaged, 
and their statistical error bars are determined. Through- 
out this work, we use the convention that one standard 
error corresponds to a one-sigma confidence interval. 

For the SiH molecule, Figure [3] shows different DMC 
force estimates defined in Eqs. (pjl- pT]) and evaluated 
at different bond lengths. For reference, this figure also 
shows the DMC energies which are simultaneously cal- 
culated with the forces. A Morse potential fitted to the 
energies is also plotted together with its gradient. The 



dotted vertical line indicates the zero-force geometry ob- 
tained from the DMC energy curve. This figure shows 
that the geometry obtained from the minimum of the po- 
tential energy curve agrees well with the zero-force geom- 
etry obtained from the pure DMC total force estimates. 



Since we investigate molecules containing H atoms 
and a heavier atom, one could obtain equilibrium bond 
lengths and vibrational frequencies from the forces on 
the H atoms alone. However, to test the force estimators 
on heavier atoms directly, we report bond lengths and 
frequencies obtained using both the zero-force condition 
on the H atom and on the heavier atoms. For SiH4, the 
estimates of the forces on the Si atom should be zero by 
symmetry and this condition is satisfied within a stan- 
dard error of less than 0.001 a. u. Also, the symmetries 
of the H2 and SiH4 molecules imply that the force on 
each H atom should have the same magnitude. Since 
we found this to hold within statistical errors, we aver- 
age symmetry-related components to further reduce the 
statistical error bar. 
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FIG. 1: (Color online) Investigation of the timestep behaviour 
of DMC forces (a.u.) and DMC energies (Ha) for the SiH 
molecule when using the small basis set. Upper graph: the 
pure DMC forces FpS*e and FpJJ as a function of the DMC 
timestep, circles indicate forces on the Si atom and crosses 
indicate forces on the H atom. A quadratic form is fitted 
to these forces. The standard errors of these forces are indi- 
cated unless they are smaller than the line width of the fitted 
functions. Middle graph: the same information as the upper 
graph for mixed DMC calculations. Lower graph: the DMC 
energy as a function of timestep for comparison, a cubic form 
is fitted to the energies. 



IV. RESULTS 

A. Definitions 

We define the errors in quantities such as bond lengths 
and frequencies as the differences between the values ob- 
tained from the forces and from the energies. 



^ymcthod ~ 2/mctliod V 



(19) 



where y is the bond length a or vibrational frequency 
CO, and the "method" can be mixed or pure DMC. E 
indicates the DMC energy and F stands for the differ- 
ent force estimators: HFT (the HFT estimator), HFT+P 
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FIG. 2: Dependence of the future-walking pure DMC HFT 
force estimates (a.u.) on the future-walking projection time 
(a.u.). Upper right graph: the HFT force on the Si atom for 
the SiH molecule as calculated within the mixed DMC and 
future-walking pure DMC method for different bond lengths 
using the small basis set. A Morse potential is fitted to these 
forces. Upper left graph: future-walking pure DMC forces cal- 
culated for different bond lengths plotted against the future- 
walking projection time. The mixed DMC force corresponds 
to zero future- walking projection time. Lower graphs: the 
same information as the upper graphs for the large basis set. 
When using the small or large basis set, we see that the future- 
walking projection time of 10 a.u. is long enough to obtain 
accurate pure DMC estimates. 



(with the pseudopotcntial Pulay estimator) and i^tot, in 
the mixed DMC method as defined in Eqs. (|6|)- (fTT|) . and 
in the pure DMC method as defined in Eqs. (fT2l) - (fT7l) . 

In mixed DMC calculations, the error Ay may arise 
from the Reynolds' approximation of Eq. (jlip and from 
replacing all total derivatives with partial derivatives. 
Both approximations introduce an error of first order in 
(v]/rp — <j)). In pure DMC calculations, the error Ay may 
arise from the approximate nodal term in Eq. p7p and 
from replacing the total with partial derivatives. The 
latter approximation introduces an error of first order in 
(vjTrp — (())^ whereas the approximate form of the nodal 
term introduces an error of second order. 



B. Bond lengths 

Table Hill presents equilibrium bond lengths calculated 
with the mixed and pure DMC total force estimators for 



-Acxpt 



'^mix \^) Q-purcj-^-j '^mix j^j 



^purc \ 



Q-nnrc \'^) 



Ha 0.741 0.7416(2) 0.74090(3) 0.7411(2) 

SiH 1.520 1.5195(7) 1.5123(1) 1.5179(8) 1.5107(1) 1.5173(9) 

GeH 1.589 1.6012(8) 1.5913(1) 1.5993(9) 1.5901(1) 1.5992(10) 

SiH4 1.480 1.4740(4) 1.46970(9) 1.4728(10) 

TABLE III: Equilibrium bond lengths (A) for four molecules calculated from mixed and pure DMC forces on each atom. The 
first column gives the experimental bond lengths from Ref.-^, the second column shows values obtained from the DMC energy 
curves, and the other columns give bond lengths obtained from the zero-force condition for the mixed and pure DMC total 
force estimators as defined in Eqs. (|6} and (|12p . We specify bond lengths obtained from forces on the H atoms by (1) and on 
the non-H atoms by (2). 



Basis 



Aa 



HFT+P/ 



1) A<-'(1) 



AaZ.W 



Aa«,^T+P(l) 



Aag'^-Ul) 



Ha large 0.7416(2) 


-0.0057(2) 


-0.0057(2) 


-0.0007(2) 


-0.0002(2) 


-0.0003(2) 


-0.0004(3) 


SiH large 1.5195(7) 


-0.0098(8) 


-0.0091(8) 


-0.0072(7) 


-0.0006(10) 


0.0007(10) 


-0.0016(11) 


GeH large 1.6012(8) 


-0.0138(9) 


-0.0141(9) 


-0.0099(8) 


-0.0012(11) 


-0.0018(11) 


-0.0019(12) 


SiH4 large 1.4740(4) 


-0.0055(6) 


-0.0056(6) 


-0.0043(4) 


-0.0005(8) 


-0.0007(8) 


-0.0012(11) 



SiH smaU 1.5242(6) 0.0006(7) 0.0023(7) -0.0091(6) 0.0052(8) 0.0086(8) 0.0004(9) 

GeH smaU 1.5991(7) -0.0059(7) -0.0058(7) -0.0060(7) 0.0024(10) 0.0028(10) -0.0010(11) 



Basis a 
SiH large 1.5195(7) 
GeH large 1.6012(8) 



Aa^U{2) Aa^f,^+P(2) Aaf^-HS) AaZ^^{2) AaZ^^+^ (2) Aagi;°U2) 



-0.0283(8) -0.0279(8) -0.0088(8) 
-0.0208(10) -0.0200(10) -0.0111(8) 



-0.0032(10) -0.0023(10) -0.0022(11) 
-0.0044(12) -0.0026(12) -0.0020(13) 



SiH smaU 1.5242(6) -0.0142(7) -0.0207(7) -0.0101(7) 
GeH smaU 1.5991(7) -0.0220(8) -0.0206(8) -0.0085(7) 



0.0233(8) 0.0090(8) 0.0016(9) 

-0.0027(10) 0.0001(10) -0.0023(12) 



TABLE IV: Differences Aa in the equilibrium bond lengths (A) derived from the DMC forces and from the DMC energies. Aa 
is defined by Eq. (|19|l . The first column specifies the basis set used, the second column gives bond lengths obtained from the 
DMC energy curves, and the other columns give bond lengths obtained from the zero-force condition for the three different 
estimators, F^^'^ , F^ft-i-p^ ^^^ ptot ^ ^^ ^-^^ mixed and pure DMC methods, defined in Eqs. ©-JlTl). We specify bond lengths 
obtained from forces on the H atoms by (1) and on the non-H atoms by (2). 



four molecules using the FPLA. Table IIVI gives further 
details of the different contributions to the total forces 
and presents the difference, Aa, of the bond lengths de- 
rived from the forces and from the energies, as defined in 
Eq. miD. 

We begin by discussing pure DMC results with the 
large basis set. The HFT estimates Fi|j^J from forces on 
the H atoms give very accurate bond lengths (upper part 
of Table HV)) . As shown in the third last column of Ta- 
ble [IVl the errors Aap^j.^(l) in the bond lengths derived 
from the forces on the H atoms are not much larger than 
one standard error of 0.001 A. The bond lengths from the 
HFT forces calculated on the non-H atoms (lower part of 
Table [IV| are not as accurate. Adding the pseudopoten- 
tial Pulay term to the HFT forces, Fi||^J+^, has a very 
small effect on the bond lengths derived from the forces 
on the H atoms, as shown in the penultimate column of 
Table IIVI For forces acting on the non-H atoms, how- 
ever, adding the pseudopotential Pulay term improves 
the bond lengths slightly. The nodal terms are very small 
and do not significantly change the bond lengths obtained 
in the large basis set pure DMC calculations. The error 
bars on Aa^^j!-^ are dominated by the contribution from 



the HFT force, so that including the pseudopotential Pu- 
lay and nodal terms does not increase the noise much. 

Table IIVI shows that both the pseudopotential Pulay 
and nodal terms become more important for SiH and 
GeH when using the small basis set. The bond lengths 
from the HFT forces on both the H and Si atoms are 
significantly worse than for the large basis set. However, 
when the pseudopotential Pulay and nodal terms are in- 
cluded the bond lengths are not significantly worse than 
for the large basis set. 

When comparing the mixed and pure DMC total force 
estimates, we find from Table Hill that the statistical er- 
rors in all bond lengths obtained from the mixed DMC 
forces arc about a factor ten smaller. This is because the 
pure DMC estimator used here does not satisfy a zero- 
variance condition. The absolute deviation in all bond 
lengths derived from the mixed DMC total forces and 
from the energies is on average 0.0076(2) A. In a sim- 
ilar comparison, the pure DMC total forces give bond 
lengths with a much smaller absolute average deviation 
of 0.0015(4) A. Although adding the Pulay terms to the 
mixed DMC HFT force may improve the bond lengths 
by up to 17 standard errors in our results, all pure DMC 
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FIG. 3: (Color online) Upper graph: three different estimates 
of the pure DMC forces, F^JJ , F™J+^ and F^e, for the Si 
and H atoms of the SiH molecule evaluated at different bond 
lengths. These estimates are defined in Eqs. (|12|l - (|17|l and are 
calculated with the small basis set. The Morse potential is 
fitted to the forces. Middle graph: the same information as 
for the upper graph with the three estimates, F™^ , F™^'^^ 
and -F^°ix, calculated within the mixed DMC method. These 
estimates are defined in Eqs. (p)l- (|ll|l . Lower graph: the DMC 
energy at different bond lengths, a Morse potential fitted to 
the DMC energies and the derivative of the Morse potential 
corresponding to the forces on the Si and H atoms. To guide 
the reader, the dotted vertical line corresponds to the zero- 
force geometry obtained from the fitted DMC energy curve. 



forces (with and without Pulay terms) generally give 
more accurate bond lengths than the best mixed DMC 
total force estimates. This difference in accuracy can be 
understood by recalling that the error introduced in the 
mixed DMC force estimates is of first order in (^t — *&) 
whereas the error in the pure DMC force estimates is 
only of second order. The additional first order error 
from replacing the total derivatives by partial derivatives 
appears to be small. 



H2 4401 4420(16) 4441(2) 

SiH 2042 2049(14) 2075(2) 

GeH 1908(35) 1907(14) 1944(1) 

SiH4 2187 2288(11) 2324(2) 



4403(15) 

2041(11) 2080(2) 2052(12) 

1909(12) 1949(2) 1904(13) 

2299(29) 



TABLE V: Harmonic vibrational frequencies (cm~^) calcu- 
lated within the mixed and pure DMC methods. The abbre- 
viations are analogous to those in Table Hill 



The differences between the DMC bond lengths (from 
either the DMC energy or the pure DMC forces) and 
experimental data in Table IIIII are somewhat larger 
than the difference between the bond lengths from the 
DMC energies and forces. This difference is largest for 
GeH (0.010(1)-0.012(1)A), followed by SiH4 (0.006(0)- 
0.007(1) A) and SiH (0.001(1)-0.003(1) A) and is negligi- 
ble for H2. These bond length deviations from experi- 
ment must largely arise from a combination of the fixed- 
node approximation, the FPLA scheme, which slightly 
alters the pure DMC ground-state distribution, and the 
pseudopotentials. The fixed-node approximation could 
be improved by using more accurate trial wavefunctions. 
It is more difficult to develop better pseudopotentials, 
although including core-polarization potentials on the Si 
and Ge atoms may also improve the result o'^^'^^i'^^ . 



C. Vibrational Frequencies 

Table |V| presents harmonic vibrational frequencies cal- 
culated with the mixed and pure DMC total force estima- 
tors. Table IVTl gives details of the different contributions 
to the total forces and presents the differences, Acj, of the 
frequencies derived from the forces and from the energies, 
as defined in Eq. P^ . 

Table IVTl shows that for all pure DMC force estimators 
the difference in the vibrational frequencies derived from 
the forces and the energies is comparable to or less than 
one standard error of about 20cm~^ with the exception 
of SiH. The effect of the pseudopotential Pulay and nodal 
terms on the pure DMC frequencies is small. This may 
also be seen qualitatively in Figure [31 where adding the 
Pulay terms mostly shifts the forces at different bond 
lengths by similar amounts. 

As in the discussion of bond lengths, we find that the 
mixed DMC total forces give less accurate frequencies 
than the pure DMC forces. The absolute difference be- 
tween frequencies derived from the mixed DMC total 
forces and from the energies is on average 43(4) cm~^ 
compared to 16(6) cm"-'^ when the frequencies are derived 
from pure DMC total forces. Although adding the Pulay 
terms to the mixed DMC HFT force may improve the 
frequencies by up to ten standard errors in our results, 
all pure DMC force estimates (with and without Pulay 
terms) still give more accurate frequencies than the best 



Basis 



auj: 



A^S-x^l) Ac.p",'-J(l) Aa.gfJ+P(l) A<i;i;(l) 



H2 large 4420(16) 128(18) 128(18) 

SiH large 2049(14) 68(15) 70(15) 

GeH large 1907(14) 87(15) 89(15) 

SiH4 large 2288(11) 71(16) 74(15) 



21(16) 


-10(19) 


-10(19) 


-16(22) 


26(13) 


-10(16) 


-7(16) 


-7(18) 


37(14) 


1(16) 


4(16) 


2(18) 


36(11) 


25(23) 


30(23) 


11(31) 



SiH small 1992(12) 71(13) 
GeH small 1905(10) 74(11) 


72(13) 
74(11) 


81(12) 
29(10) 


16(14) 
-25(13) 


19(15) 
-25(13) 


46(16) 
0(15) 


Basis o)^ Aa;Sfx^(2) 


Aa-:+-(2) 


AaS°.'(2) 


Aaflli2) 


Aa«r-+P(2) 


Aa^-U2) 


SiH large 2049(14) 103(16) 
GeH large 1907(14) 73(16) 


100(17) 
68(16) 


31(14) 
41(14) 


14(16) 
2(17) 


1(17) 
-10(17) 


3(18) 
-4(19) 


SiH small 1992(12) 104(13) 
GeH small 1905(10) 86(12) 


129(13) 
85(12) 


92(12) 
38(10) 


-23(17) 
23(14) 


27(15) 
17(14) 


53(16) 
20(16) 



TABLE VI: Difference Alo in the harmonic vibrational frequencies (cm ^) derived from the DMC forces and energies. The 
difference Au is defined by Eq. (|19|l . The abbreviations are analogous to those in Table HVl 



mixed DMC total force estimates. 



D. Comparison of the FPL A and SPLA schemes 

The FPLA and SPLA schemes are compared when cal- 
culating forces for the SiH and GeH molecules. Since 
the schemes generate different ground-state wavefunc- 
tions, expectation values may also differ between the two 
schemes. Additionally, we also use slightly different ap- 
proximations in the force estimators under the FPLA and 
SPLA schemes. The calculated bond lengths and vibra- 
tional frequencies are presented in Tables IVIII and IVIIII 
When bond lengths and frequencies are derived from any 
of the three pure DMC force estimates, the results be- 
tween the two localization schemes agree within or close 
to one standard error of about 0.0015 A for the bond 
lengths and about 20 cm^^ for the frequencies. When the 
results are instead derived from the mixed DMC forces, 
we find that the difference can be as large as three stan- 
dard errors. 



V. CONCLUSION 

We have presented exact expressions for the total 
atomic forces within mixed and pure diffusion Monte 
Carlo (DMC) calculations using nonlocal pseudopoten- 
tials and reported approximations for estimating them. 
These expressions include the pseudopotential Pulay 
termi^ and the Pulay nodal tcrmi. 

We obtained harmonic vibrational frequencies and 
equilibrium bond lengths from the calculated forces for 
four small molecules. The calculations were performed 
with single-determinant Slater- Jastrow trial wavefunc- 
tions using the mixed DMC and future-walking^ pure 
DMC methods. In the pure DMC calculations we found 
that the contributions to the force from the Pulay nodal 



term and the pseudopotential Pulay term were compa- 
rable to or less than the statistical error in the total 
force, when the trial wavcfunction and the nodal sur- 
face were sufficiently accurate. In these cases, neglecting 
the nodal and pseudopotential Pulay terms could have 
been justified. However, when the trial wavefunctions 
were less accurate, both Pulay terms became important 
and including them significantly improved the pure DMC 
forces. All bond lengths and vibrational frequencies de- 
rived from the pure DMC total forces agreed with those 
obtained from the DMC energies within or close to one 
standard error. This showed that the error from replac- 
ing total with partial derivatives in the pure DMC force 
estimators is very small, and that the additional error 
from approximating the pure DMC nodal term is well 
behaved. 

The deviations of the bond lengths and frequencies ob- 
tained from the mixed DMC total forces and from the en- 
ergies were generally much larger than in the pure DMC 
calculations. This was explained by the less accurate ap- 
proximations in the mixed DMC force estimates which 
introduce errors of first order in (^^t — $)■ For a speci- 
fied quality of trial wavcfunction we therefore concluded 
that pure DMC forces were more accurate than mixed 
DMC ones. We also investigated both the FPLA and 
SPLA schemes for treating the nonlocal pseudopotential 
operator and found very similar results. 

A brief review of previous attempts to calculate forces 
within the DMC method and a discussion of the perfor- 
mance of various quantum chemistry methods in estimat- 
ing bond lengths and vibrational frequencies for several 
molecules was presented in RefJ^. The deviation be- 
tween our results and experimental data is comparable to 
or better than results obtained by other accurate quan- 
tum chemistry methods, and is generally considerably 
better than in density functional methods. Our work has 
demonstrated that accurate atomic forces can be calcu- 
lated with pseudopotcntials and the DMC method. 
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PLA 



a^ Aa^L^(l) Aa::^f,^+^(1) Aa^r^JX) Aa^^^il) 



Aa 



HFT+P 



(1) Aafi-Ul) 



SiH FPLA 1.5195(7) -0.0098(8) -0.0091(8) 

SPLA 1.5210(8) -0.0113(9) -0.0110(9) 

GeH FPLA 1.6012(8) -0.0138(9) -0.0141(9) 

SPLA 1.5995(9) -0.0129(10) -0.0128(10) 



-0.0072(7) -0.0006(10) 
-0.0086(8) -0.0024(10) 
-0.0099(8) -0.0012(11) 
-0.0073(9) -0.0004(11) 



0.0007(10) -0.0016(11) 

-0.0019(10) -0.0034(11) 

-0.0018(11) -0.0019(12) 

-0.0004(11) -0.0039(13) 



PLA a^ AaSfJ(2) Aa^^:^+^ (2) Aa^^{2) AaZl{2) Ao 



'{2) Aa^^°U2) 



SiH FPLA 1.5195(7) -0.0283(8) -0.0279(8) 

SPLA 1.5210(8) -0.0303(9) -0.0289(9) 

GeH FPLA 1.6012(8) -0.0208(10) -0.0200(10) 

SPLA 1.5995(9) -0.0157(10) -0.0148(10) 



-0.0088(8) -0.0032(10) 
-0.0098(8) -0.0055(10) 
-0.0111(8) -0.0044(12) 
-0.0084(9) -0.0060(12) 



-0.0023(10) -0.0022(11) 

-0.0026(10) -0.0034(12) 

-0.0026(12) -0.0020(13) 

-0.0042(12) -0.0033(15) 



TABLE VII: Difference Aa in the equilibrium bond lengths (A) derived from the DMC forces and energies for SiH and GeH. 
The FPLA and SPLA localization approximations are used as indicated in the first column. All other abbreviations are the 
same as in Table HVl 



PLA 



^" Ac.g!7(l) Ac.;^|:^'+^(1) A^;Sx (1) 



A<^J(1) Ac.fr^+^{1) A<SU1) 



SiH FPLA 2049(14) 68(15) 70(15) 26(13) -10(16) -7(16) -7(18) 

SPLA 2050(12) 73(14) 73(13) 24(12) 13(15) 11(15) 6(18) 

GeH FPLA 1907(14) 87(15) 89(15) 37(14) 1(16) 4(16) 2(18) 

SPLA 1915(13) 33(14) 33(14) 20(13) -19(15) -19(15) 20(18) 



PLA 



A<f7(2) A<^^^+^(2) Aa^,°;(2) Aa^„^;^'(2) Aag,^;^+"(2) Aa-^{2) 



SiH FPLA 2049(14) 103(16) 100(17) 31(14) 14(16) 1(17) 3(18) 

SPLA 2050(12) 95(16) 89(15) 28(12) 19(15) 6(16) 2(19) 

GeH FPLA 1907(14) 73(16) 68(16) 41(14) 2(17) -10(17) -4(19) 

SPLA 1915(13) 25(15) 25(15) 24(13) 3(17) -3(16) 2(20) 

TABLE VIII: Difference Alu in the harmonic vibrational frequencies (cm^^) derived from the DMC forces and from the DMC 
energies for SiH and GeH. The FPLA and SPLA localization approximations are used as indicated in the first column. All 
other abbreviations are the same as in Table HVl 
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